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1  Introduction 


In  this  paper  we  introduce  a  new  spectral  technique  for  the  numerical  solution  of  nonlinear 
hyperbolic  equations.  This  technique,  as  almost  every  modern  finite  difference  scheme  for 
shock  computations,  is  based  on  the  cell  averaged  form  of  the  equations.  This  is  essential  for 
finite  difference  shock  capturing  techniques  and  it  is  our  experience  that  it  plays  an  essential 
role  in  a  successful  spectral  simulation  of  problems  that  involve  shock  waves  [1]. 

Consider  the  nonlinear  hyperbolic  equation 


{ 


ut  +  fx(u)  =  o,  x  e  [-1,1] 

U(x,0)  =  Uo(x), 


(1.1) 


with  appropriate  boundary  conditions. 

The  cell  averaged  form  of  (1.1)  is  obtained  by  integrating  (1.1)  between  any  two  points 
— 1  <  a  <  b  <  1  to  get 


-  irb  f  u{x)  **  -  - F{ 'u^-  <12> 

Let  u(x,i)  be  an  approximation  to  U(x,t)  at  time  t.  Following  Harten  [4]  we  express  the 
approximation  to  U(x,t  +  r)  by 


u(x,t  +  t)  =  AE(t)  %{• ;  ■&), 


(1.3) 


where  A  is  the  cell  averaging  operator  and  E(t)  is  the  exact  time  evolution  operator  corre¬ 
sponding  to  (1.1).  Throughout  the  paper  we  will  not  distinguish  between  7Z(-;u)  and  7iu 
The  operator  7t(-;  it)  is  of  extreme  importance,  it  represents  the  way  we  reconstruct  u  from  its 
given  cell  average  values  Uj_  i  =  S*j -l  u(x)  dx,  where  {x^}^  are  the  grid  points.  For 

finite  difference  schemes  u  is  a  piecewise  polynomial  of  low  degree,  so  that  the  reconstruction 
itself  is  simple.  It  becomes  complicated  only  if  one  imposes  also  the  requirement  that  the 
reconstruction  should  be  essentially  nonoscillatory.  In  [1]  we  have  presented  an  essentially 


non-oscillatory  Fourier  method  based  on  the  cell  averaging  formulation  (1.2).  In  that  case  codea 


jd/or 

Special 
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Jlat 


□  □ 


the  transformations  between  the  cell  averages  and  the  point  values  are  simple  and  can  be 
carried  out  efficiently  by  the  Fast  Fourier  Transformations  (FFT).  This  can  be  attributed  to 
the  fact  that  the  boundary  conditions  are  periodic  and  that  the  cell  average  of  a  trigonomet¬ 
ric  function  is  proportioned  to  the  function  itself.  However,  for  Chebyshev  methods,  the  cell 
averaging  operation  (denoted  by  the  operator  A)  is  not  simple  nor  is  %(■ ;  it).  As  a  matter  of 
fact  not  only  the  formulation  but  edso  the  implementation  of  A  eind  'R,  is  not  straightforward. 

In  this  paper  we  formulate  the  cell  averaging  technique  for  the  Chebyshev  method.  We 
will  discuss  its  stability  for  linear  problems  and  show  an  example  of  its  applicability  to 
nonlinear  systems  of  equations  by  simulating  the  problem  of  shock-density  wave  interaction. 
The  cell  averaging  formulation  is  an  essential  part  of  the  numerical  code. 

The  outline  of  the  paper  is  as  follows:  In  Section  2  we  show  how  to  reconstruct  efficiently 
point  values  of  a  polynomial  from  its  cell  averages  and  vice  versa.In  Section  3  we  introduce 
the  new  numerical  technique  and  show  its  stability  for  linear  problems.  Section  4  is  devoted 
to  numerical  results  obtained  by  using  the  new  method. 

2  Cell  Averages  and  Point  values 

In  this  section  we  will  discuss  the  cell  averaging  operator  A  and  the  reconstruction  operator 
%  in  the  context  of  the  Chebyshev  methods.  In  these  methods  the  approximations  are 
taken  from  the  space  of  polynomials  of  degree  N.  It  is  therefore  clear  that  A  and  H,  when 
restricted  to  the  polynomial  space,  can  be  expressed  as  matrices  Ajy  and  Rjy.  We  will  give 
the  explicit  formulas  for  these  matrices.  We  start  by  discussing  the  operator  A. 

Assume  that  f{x)  is  in  CT[— 1, 1],  r  >  0.  Let  Xj  =  cos  (^),  0  <  ;  <  N  be  the  Chebyshev- 
Lobatto  points  in  [-1,1].  For  later  use,  define  x^_i  =  cosf  -),  1  <  j  <  A.  The  cell 
averaged  function  f(x)  of  f(x)  is  defined  as 

1  rfj(x) 

/(I)  “  Af  =:  A.(.) /(I)  ix  fOT  -1  < 1  <  ‘'  <21> 


2 


where 


,  ,  x  ,  _!  A0 

/ii(x)  =  cos  (cos  i - — 

,  ,  ,  .  A 0. 

/i2(s)  =  cos  (cos  x  +  — ), 

4 

.  .  7T 

A*  =  T 


(2.2) 


The  reason  for  the  definition  in  (2.1)  is  that 

f(xj-0  =  ~ — " — T*  /  '  f(x)dx>  for  1  <  j  <  N.  (2.3) 

Xj—l  Xj  Jxj 

As  stated  before,  we  are  interested  in  A  operating  in  a  polynomial  space.  In  Lemma  1, 
we  show  that  the  result  of  A  on  a  polynomial  is  still  a  polynomial  of  the  same  degree. 


Lemma  1  LetUk(x)  =  j^T^^x),  k  >  0,  be  the  second  kind  of  the  Chebyshev  polynomials. 
Then 

Uk(x)  =  AUk  =  akUk(x),  (2.4) 

where 


sin  (k  +  l)x 
( k  +  1)  sin  2 


ak  = 


as  • 


(2.5) 


Proof:  Substituting  Uk(x)  in  the  right  hand  side  of  (2.1)  and  making  the  transformation 
x  =  cos  6,  0  <  0  <  7r,  we  have 


Uk(x)  = 


COS  (cos-1  X  +  A£)  —  COS  (cos-1  X  -  ^p)  Jco»{ CO*-1 1-^2) 


3-1x_M)/, 


co§  (co*-1  x+£~) 


C/jt(x)  dx. 


Since  T*(x)  =  cos  (k6),  then  Uk(x)  =  and  therefore 


1  re+¥ 

t4(x)  =  -zr, - tt—. — a T~- — n  /  sin(A;  +  l)0d0 

v  '  2(A:  +  l)sin  A£  sin0  v  ' 


1 


2(k  +  l)sin  A£  sin0 
sin(fc  +  l)--9 


,  1%cos(A:  +  1)0 
1  '  k+l  |9-¥ 


( k  +  l)sin 


1)^T  /  1  sin(k  +  1)0 \ 

in  ~  \  k  +  1  sin  0  J  ’ 
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i.e. 


Uk{  x)  =  akUk(  x). 


Q.E.D. 

i,Prom  Lemma  1  one  gets 

Corollary  1  The  cell  averaged  function  of  any  polynomial  of  order  N  is  a  polynomial  of  the 
same  order. 


Proof:  any  polynomial  of  degree  N  has  the  expression 

f(x)  =  akUk{x ), 

k=o 


where  ak  are  constants. 
Therefore,  by  Lemma  1 


(2.6) 


/(*)  =  E  cr»a*£4(1)'  (2.7) 

fc=0 

Q.E.D. 

Thus,  Lemma  1  gives  explicitly  the  eigenvalues  <7*  of  the  matrix  Ajv  (the  restriction  of 
A  to  polynomials  of  degree  N).  These  eigenvalues  are  uniformly  bounded  from  above  and 
below,  in  fact 

-  <  C7fc  <  0  <k<N.  (2.8) 

7T  2 

If  /( i)  is  a  polynomial  of  degree  N,  then  it  is  uniquely  determined  by  its  values  f(xj),  j  = 
0,  •  •  • ,  N.  So  theoretically  f{Xj_  i ),  j  =  1,  •  •  • ,  N  can  be  determined.  Therefore  the  transfor¬ 
mation  from  f(xj),  j  =  0,  ■  •  • ,  N  to  f(xj_  i),  j  =  1,  ■  •  • ,  N  is  well  defined  and  we  only  need 
to  address  the  issue  of  its  efficient  implementation. 

In  general  it  is  known  that 

/(*)  =  E  onH(»).  (2.9) 

fc= 0 
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where 


“  jy-^  f(.Xi)Tk{Zj), 

cq  =  cn  =  2,  ck  =  1  if  k  ^  0  or  N.  (2.10) 


Alternatively 

/(®)  =  £  *>kUk(x),  (2.11) 

k=o 


bk  =  y,  *  =  N  -  1  or  N, 

h  =  ^( ckak  -  ojfc+a),  0  <  k  <  N  -  2  (2.12) 


and  therefore  by  corollary  1 

Kx)  =  J2akbkUk(x).  (2.13) 

k=0 

Now  {/(Sj-^J-yLi  are  obtained  by  substituting  Xj_ ^  in  (2.13)  (this  can  be  carried  out 
using  the  FFT). 


We  note  that  equations  (2.9)  -  (2.13)  describe  how  to  get  the  vector  f(x^_  i ),  j  —  1,  •  •  • ,  N 
from  the  vector  f(xi),  i  =  0,  •  •  • ,  N.  Denote  by  An  the  N  x  (N  +  1)  matrix  defined  by  this 
transformation.  We  note  that  An  can  be  written  explicitly.  In  fact  the  polynomial  /(x)  has 
an  unique  representation  as 

/(*)  =  £  /(*;)&(*)>  (2-14) 

3=0 

where 

9'(  >  a,JV’(x  -  x,)  ctN  to  ft 

with  cfc  defined  in  (2.10)  . 

It  follows  upon  substituting  (2.14)  in  (2.1)  and  using  the  fact  that  T0(x)  =  Uq{x),  T\(x )  = 
2<j\U\(x)  and  Tk(x )  =  \{akUk{x)  -  crfc_at4_ a(x))  for  k  >  2, 

/(*)  =  E/(*<)&(X)»  (2.15) 

l=o 
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where 


3t(x)  =  \  +  2ctiTi(x/)[/i(x)  + 


2  ^  rfc(a£)(gfcyfc(a)  -  <rk-3Uk-a(x)) 
ctN  ~  2cfc 


with  <7*  defined  in  (2.5). 
Setting  x  —  Xj_  i  in  (2.16) 


(2.16) 


(AN),<  =  1),  0<l<N,l<j<N.  (2.17) 

Thus  we  have  outlined  two  procedures  to  get  /(x3_ i )  from  f(xj),  one  uses  the  FFT  and 
another  uses  matrix  vector  multiplications. 

We  are  now  ready  to  discuss  the  reconstruction  operator  7£(-;/).  Note  that  in  the 
beginning  of  the  solution  process  (1.3),  we  only  have  the  values  /3_i,  j  =  1,  •  •  • ,  N,  thus  we 
need  another  piece  of  information  in  order  to  define  uniquely  the  AT-th  degree  polynomial 
f(x).  This  piece  of  information  is  provided  by  the  boundary  condition.  For  simplicity,  we 
assume  that  the  boundary  condition  is  of  the  form 


/( 1)  =  /(*o)  =  fo- 


(2.18) 


The  reconstruction  is  done  in  two  steps.  Define  first  a  (N  —  l)-th  polynomial  ftf-i(x) 
which  collocates  /(»)  at  {x3_i}f=1  ,  i.e.  /*-i(*,_*)  =  /(*,-*).  1  <  i  <  ^  is  ^adily 

verified  that 

fN-i(x)  =  £  ckTk(x),  (2-19) 

k=0 

where 

ck  =  -^rry)/(x,  i)Tfc(x  _i),  cfc  is  same  as  in  (2.10).  (2.20) 

5  3 

Alternatively  we  have 

/».,(») = r  ww.  <2-21) 

fc=0 

where 


bk  =^,  k  =  N  —  2,  AT  —  1, 

ifc  =  2^kCk  ~  c*+2)>  0  <  k  <  N  -  3. 


(2.22) 


Now,  by  Corollary  1 


h-i(x)  =  £  —Uk(x)  (2.23) 

k—0  a* 

is  a  polynomial  of  degree  N  —  1  such  that  A/n-i  =  /jv-i- 

Generally,  fN-i(x)  does  not  satisfy  the  boundary  condition  (2.18).  There  are  two  ways 
to  modify  /n_i(x)  so  that  the  boundary  condition  (2.18)  is  satisfied.  In  the  first  way  we  can 
add  to  fa. i(x)  an  IV-th  degree  polynomial  Q(x )  such  that 

$(*;-*)  =  0,  1  <j<N,  (2.24) 

and 

In- i(a)  +  Q{  1)  =  /o-  (2.25) 

It  can  be  verified  that  in  order  to  satisfy  (2.24) 

Q(x)  =  c((l-x,)T'„(x))'.  (2.26) 

Let  /(x)  be  the  sum  of  /n-i(s)  and  Q{x), 

/(«)  =  /*-,(*) +  c((l-*1)r„(x))'  (2.27) 

=  fN-1  (x)  -  c[xT'n(x)  +  N2Tn(x)]. 

The  last  equality  follows  from  the  Chebyshev  equation,  and  the  constant  c  is  now  deter¬ 
mined  by  the  condition  f(x 0)  =  /0,  i.e. 

c=-2^(/o-/*-i(1))-  (2.28) 

Finally  given  f{Xj_  i ),  j  =  1,  •  •  • ,  N  and  fo  we  can  get 

f(xi)  =  /tf-i(xi)  -  c  \xiT'N{xi)  +  JV3rw(xi)]  ,  i  =  0,  ■  •  • ,  N  (2.29) 

where  can  be  evaluated  from  (2.23)  by  using  the  FFT. 

Note  that  in  this  procedure  we  change  the  values  of  at  all  the  grid  points. 


n 

i 


Denote  by  Rn  the  ( N  +  1)  x  ( N  +  1)  matrix  transforming  f0  and  /( 1  <  j  <  N 
to  f{xj),  0  <  j  <  N,  i.e. 


(/M, •  •  • , I{xn))J  =  Rn  (/(so), /(si), •  •  • , /(Stf.i ))T  . 

As  before  we  can  write  Rn  explicitly.  Equation  (2.19)  can  be  rewritten  as 

In- i(s)  =  S/(s,-^)^i(s), 

3=1 

where  hj  (sfc_i)  =  Sjk  and  hj(x)  are  polynomials  of  degree  N  —  1;  explicitly 

**(*> =  (-tf*1* (a  -  5>^)^rT- 

J- 5 

It  can  be  shown  that  /ij(x)  is  the  cell  averaged  polynomial  of 


hj(x)  =  £  ^kUk(x) 
k=o 


with  A*,  defined  as 


(2.30) 


(2.31) 


(2.32) 


(2.33) 


1 


Xk  =  N-2,N-l, 


At  = 


Nak 

1 


(T^.i)  -  r*+a(*j-_i)),  if  o  <  fc  <  AT  -  3  .  (2.34) 


As  a  result  of  (2.31)  and  (2.33)  polynomial  fN~\(x)  takes  the  form 


JV-1 


/w-i(s)  =  £  /(**-*)*»(*)» 


J=1 


and  by  (2.27) 


/(*)  =  /»-.(*)  ~(Ar  /»-.(!))  ((1  -  *’»(»))'• 


(2.35) 


Using  (2.31)  we  get 

/(*)  =  -  ftp  (/» ((i  -  *’»(*))'  (2.36) 

=  e/(»h>  M*) + wr  ((!  -  x’ra*))'!  -  (d  -  *’»(*))' 
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Substituting  x  =  X{  into  (2.36)  gives 


'  8 

(Rn)„  =  I 


if  i  =  0, 
if  j  =  0, 
l  <  i,  j  <  N. 


(2.37) 


To  summarize  we  state 


Lemma  2  Let  f(x)  =  7£(- ;/)  be  the  N-th  polynomial  defined  in  (2.27)  then 


i)  =  f(xj_  i)  /or  1  <  j  < 


N. 


(2.38) 


A  different  way  to  modify  /jv_i(a:)  in  (2.23),  in  order  to  satisfy  the  boundary  condition 
(2.18),  is  to  add  to  it  a  polynomial  Qi(x )  of  degree  N  such  that  the  point  values  /v-i(®i) 
remain  unchanged  and  the  new  polynomial  satisfies  the  boundary  condition.  Thus  instead 
of  (2.27)  we  define 

/(*)  =  /»-.(*)  +  (A  -  /*-.(l))(1  +2^?'(X).  (2  39) 

Computationally  (2.39)  is  simpler  than  (2.23)  (2.27).  However,  note  that  in  this  case 

^  (2-40) 

which  is  in  contrast  to  (2.38).  The  matrix  corresponding  to  (2.40)  can  be  formed  similarly 
as  in  (2.37). 


3  Cell  Averaging  Chebyshev  (CAC)  Method  and 
Linear  Stability 

In  this  section  we  will  establish  the  stability  of  the  Cell  Averaging  Methods,  presented  in 
Section  2,  when  applied  to  a  first  order  scalar  hyperbolic  equation.  It  is  tempting  to  try  to 
obtain  stability  in  the  L1  norm  because  of  the  way  the  method  (1.3)  is  presented.  However 
we  will  only  give  the  stability  estimate  based  on  a  weighted  Chebyshev  norm. 
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Consider  the  initial  boundary  value  problem  of  the  scalar  hyperbolic  equation 


Ut  =  Ux,  a  6  [-1,1], 

«  U(x,0)  =  Uo(x),  (3.1) 

u{i,t)  =  o. 

The  cell  averaged  form  of  (3.1)  is 

wAx’  0  +  TT^r-TT^  PM*)’  0  -  £«*).  01  =  0.  (3.2) 

ot  h2\x)  hi(x) 

where  hi(x )  and  /^(x)  are  defined  in  (2.2). 

We  follow  the  notation  (1.3).  Let  the  N- th  degree  polynomial  u(x,i)  be  the  approxima¬ 
tion  to  (3.2).  Our  aim  is  to  find  the  error  equation  for  lZ{-\u)  defined  either  by  (2.27)  or 
(2.39).  i,From  (1.3'' 

u(x,t  +  t)  =  AETZu(x,t).  (3.3) 

Applying  the  reconstruction  operator  1Z  we  get 

7Zu(x,t  +  r)  =  TlAEHu(x,t)  (3.4) 


iFrom  equation  (3.4)  it  is  clear  that  Hu  satisfies  exactly  the  equation 


dHu  dHu  r 

-ar  =  -ar  +  rl(1 


for  the  reconstruction  procedure  (2.27)  and 


dTZu 

~W 


dl Zu 
dx 


+  Ti(l  +  x)Tff 


(3.5) 


(3.6) 


for  (2.39).  r  and  T\  are  quantities  depending  on  time  t. 

Note  that  (3.6)  is  the  error  equation  for  the  Chebyshev  collocation  method.  For  a  scalar 
linear  equation  the  CAC  method  corresponding  to  (3.6)  is  equivalent  to  the  known  collocation 
method  [2].  It  remains  to  investigate  the  stability  of  (3.5).  For  simplicity,  in  the  remainder 
of  this  section  we  denote  7Zu  by  u.  From  the  construction  of  1Z  in  section  2  we  know  that 
u(x,t)  satisfies  the  boundary  condition  in  (3.1),  i.e.  u(l,t)  =  0. 


10 


It  is  interesting  to  note  from  (3.6)  that  the  CAC  method  with  reconstruction  operator 
(2.27)  is  equivalent  -  for  constant  coefficient  case  -  to  the  collocation  method  where  the  grid 
points  are  the  zeros  of  the  polynomial 


<?(*)=  [(1 -*’»(*)]'.  (3.7) 

Note  that  by  using  standard  identities 

Q(x)  =  -  [xT;(x)  +  N2Tn(x)\  .  (3.8) 

The  term  r  in  the  right  hand  side  of  (3.5)  is  determined  by  the  boundary  condition 
u(l, t)  =  7£u(l,f).  In  fact  substituting  x  =  1  in  the  equation  (3.5)  and  noting  that  |^(1, t)  = 


0  we  get  for  r 


u*(M) 


An  alternative  expression  can  be  obtained  by  equating  the  highest  coefficient  on  both 


sides  of  (3.5),  thus  if  u  is  expanded  as 


u(x )  =  WfcTjfe(x) 

fc=0 


(3.10) 


—  1  dutf  , 

T~~N(N  +  1)  dt  1 

Before  stating  the  main  stability  result  of  this  paper,  we  state  the  following  lemma 


(3.11) 


Lemma  3  If  f(x)  —  Y^i=o  1  akTk(x),  then 


n  ^  1  ..  .  f 1  fix) 

N^T/{xj)  =  LvT^dx  +  za3N 


(3.12) 


where  Cj  is  defined  in  (2.10);  Xj  0  <  j  <  N  are  the  Chebyshev-Lobatto  points. 


Proof:  see  [3]. 

The  main  stability  result  follows  from  the  following  lemma 
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Lemma  4  Let  U(x,t )  be  the  solution  of  (3.2)  and  u(x,t)  =  TZu  be  the  approximation  to 
U(x,t )  by  the  C AC  method  (1.8)  with  the  reconstruction  operator  (2.27)  or  (2.30),  then 


dt 


7T  y,  1  1  +  Xj  2 
2 N  1  —  XjU 


tN 


2{N  +  l) 


u2N(t) 


<  0. 


(3.13) 


Proof:  We  multiply  (3.5)  by  ^j^ju(xj,t)  and  sum  from  the  points  x^  =  —  1  to  xi  to  get 

d  V  ^  1  1+Xj  2/_  7T  ^  1  l  +  x,-  ^ 

dt2N^[cjl-xju('Xj,t)  ~  N^cj  l_Xju^t>^x^ 


T2  n  r-.  1  1  +  Xj 


-  (3.14) 

iV  j=l  Cj  1  “  Xj 

We  treat  the  two  terms  in  the  right  hand  side  of  (3.14)  separately. 


1  = 

JV  Cj  i  Xj 

=  TrE  0  +  |r«J(l,  i). 

•'v  i=o  ci  1  ~  ■'v 

By  noting  the  exactness  of  the  Gauss  Lobatto  formula  one  gets 

1 = L  + ^u*(i’ () 

and  integration  by  part  yields 


=-/: 


1  l+f  u»(i,Q 


7r^(T^r*+:r’(M)' 


We  use  again  the  Gauss  Lobatto  formula  to  get 

7T 


*  l 


and  therefore 


t _  n  v'  i  n  i  xi\ u2(x}’t)  |  ^  2/-1 


We  turn  now  to  the  second  term  in  (3.14) 

II 


=  ~TN3lj  £  ^■T^~u(xLt)TN(xj) 

it  J=1  Cj  1  Xj 


(3.15) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 


(3.20) 


II=-rN>^fi^\±^x,,tyrK(xi)-TNt^l,t).  (3.21) 

Using  (3.9)  and  (3.11)  for  r  one  gets 

N  dilfj  7T  1  1  +  X;  ,  .  _«  .  .  7T 

11  =  NT\~dTN% 1)Tn{X])  ~  2NUxi-l,t)-  (3-22) 

As  ^|u(x,i)Tjv(x)  is  a  polynomial  of  degree  less  than  AN  —  1,  one  gets  using  lemma  3 

*  £  1  1  +  .V7W-,  ^  _  f1  1  +  2  uT*  ^  , _  /o  «v^ 


Toll  +Xj  .  ,  .  f1  1  +  X  UiN  , 

T7  X,  7:7— iw(xi,  <)7V(s,)  =  /  •_  —  dx  +  7ra2W 

TV  j_ q  Cj  1  Xj  J — 1  1  x  —  x 

where  a2#  is  the  2N  coefficient  in  the  expansion  of  |^u(x,  t)Ttt(x) 

It  can  be  easily  verified  that 


(3.23) 


/i&T&Kfc  =  -!*»’ 


02Af  =  —\^N, 


(3.24) 

(3.25) 


rr  W  7T  2.  . 

77  =  -’rFTT“J''flr  -  2jvu'(1,t)' 


Using  (3.19)  and  (3.26)  we  get 


t  ,  tt  ^  w  2/1 

7  +  77--’rATTI,"'-r-4]7“'(1'i)- 


(3.26) 


(3.27) 


Substituting  (3.27)  into  (3.14)  yields  (3.13).  This  completes  the  proof  of  the  lemma. 


Q.E.D. 


We  note  that  from  Lemma  3 


7T  1  1  +  Xj  2f  j\  1  y1  1  +  X  U2(x,t)  j  ^  ^ 

w  g  *>  -  2  L  —xvrr?' dx  ~  4 


(3.28) 


so  that  we  can  finally  state  that 


Theorem  1  Let  u(x,<)  =  IZu  be  the  solution  of  the  CAC  method  with  the  reconstruction 
operator  (2.27)  or  (2.  SO),  then 

1  />  1  +  x  u’(x,t)  ,  ,  2J\T  —  1  ..l  ,1  1  +  x  u2(x,0)  ,  ,  2JV  -  1 


1  /•  1  +  x  u‘  x,t  21V  -  1  2  1  /■>  1  +  x  u-  x,0  , 

2  L  l-xVl-x2<iX  +  4(¥TTj’rUK(i)  5  2  L  1-xVl-x2  *  + 


2W-1 
4(W  +  1) 
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iru^(O). 

(3.29) 


4  Numerical  Results 


In  this  section  we  apply  the  C AC  spectral  method  (1.3)  to  the  one  dimensional  gas  dynamics 
equations.  The  time  evolution  is  done  by  the  Runge-Kutta  type  method.  Each  step  of  the 
Runge-Kutta  scheme  is  done  as  follows: 


Fully  Discretized  CAC  method 
Step  1:  Reconstruction: 

given  u"  ,,  j  =  1,  •  •  • ,  N,  we  use  the  boundary  condition  and  the  matrix  Rn  to  find  the 

}~3 

point  values  uj,  j  =  1,  •  •  •  N  as  suggested  in  (2.27)  or  (2.30); 

Step  2:  Solution  in  time: 

we  update  the  values  ,  j  =  1,  •  •  • ,  N  using  the  forward  scheme, 


:~k  ~  v 


A  t 


»  xi~  i 


- -F(«?)l  .  J  =  !.•••. H, 


(4.1) 


where  A t  is  the  time  step. 

The  reconstruction  operator  Rn  yields  spectrally  accurate  point  value  approximations 
to  the  exact  solutions  if  the  exact  solutions  are  smooth,  thus  Uj_  i,  j  =  can  be 

expected  to  approximate  their  cell  averages  accurately.  However,  if  the  exact  solutions  are 
discontinuous,  the  point  values  Uj,  j  =  0,  •  •  • ,  N  obtained  by  Rn  will  be  oscillatory  as  the 
result  of  the  Gibbs  phenomenon.  In  [1]  we  proposed  a  practical  way  to  obtain  essentially 
nonoscillatory  spectral  reconstruction  to  a  discontinuous  function  from  its  oscillatory  Fourier 
approximations.  The  key  idea  there  is  to  augment  the  Fourier  space  by  adding  simple  dis¬ 
continuous  functions  whose  locations  and  magnitudes  of  discontinuities  are  approximations 
to  those  of  the  shock  waves  in  the  numerical  solutions.  In  our  computations  of  CAC  method, 
we  extend  this  idea  along  the  same  line  to  obtain  essentially  nonoscillatory  reconstructions. 
The  estimates  on  these  reconstructions  will  be  appearing  in  a  separate  work.  We  refer  the 
reader  to  ([1])  for  further  details. 

Now  consider  the  Riemann  problem  for  the  Euler  equations  for  a  polytropic  gas 
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Ut(x,0  +  4(U)  =  0, 


(4.2) 


where  U(x,  t)  and  f(U)  are  defined  as 

U  =  (p,  M,  E)T ,  — 1  <  x  <  1, 

f(U)  =  9U  +  (0,P,,P)T,  (4.3) 

where 

P  =  (7  -  1)(^  -  J/V),  W  =  W,  (4.4) 

with  7  =  1.4  and  the  initial  conditions  are  as  follows 

(plAl,Pl)  =  (3.857143,2.629369,10.33333)  when  x  <  -0.8, 

(pr,Qr,Pr)  =  (1  4-  esin57rx,0, 1)  when  x  >  —0.8,  (4.5) 

where  e  =  0.2. 

The  solutions  to  (4.2)  -  (4.5)  model  the  interaction  between  a  moving  shock  wave  and 
disturbances.  Note  that  in  the  right  state  of  the  density  a  sinuous  perturbation  of  magnitude 
e  =  0.2  is  superimposed  upon  a  constant  state.  ^From  linear  analysis  it  can  be  shown  that 
the  disturbances  will  interact  with  the  shock  wave.  A  density  wave  of  different  frequency 
and  magnitude  will  emerge  behind  the  shock  wave.  Also  the  disturbance  in  the  density 
field  will  perturb  the  velocity  and  pressure  fields  behind  the  shock  wave.  The  numerical 
solution  of  this  Riemann  problem  mandates  a  high  order  scheme  in  order  to  capture  the  fine 
structures  in  the  solutions  for  a  correct  interpretation  of  the  physical  process.  We  test  this 
problem  with  second  order  MUSCL  scheme  [5]  and  third  order  point  value  version  ENO  finite 
difference  scheme  [6]  and  the  CAC  spectral  method  proposed  in  this  paper.  Our  numerical 
results  have  shown  clearly  the  advantage  of  a  higher  order  numerical  scheme  for  problems 
with  complicated  structures. 

Figure  1-3  show  the  density  profiles  obtained  by  the  three  methods  mentioned  above. 
Figure  1  is  the  result  using  CAC  spectral  method.  The  solid  lines  are  the  solutions  taken 
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from  [6]  using  the  third  order  ENO  finite  difference  method  with  N  =  800  which  we  take  as 
the  exact  solutions  .  Figure  2  and  Figure  3  are  the  results  with  the  second  order  MUSCL 
scheme  [5]  and  the  third  order  ENO  finite  difference  scheme  respectively.  In  all  three  cases 
we  use  the  same  amount  of  mesh  points  N  =  220.  All  the  results  are  plotted  at  the  same 
time  t  =  0.3. 
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Figure  1:  The  CAC  spectral  method:  density,  N  =  220,  time  t  =  0.3.  (+)  -  numerical 
solutions,  solid  line  -  exact  solutions. 


18 


Figure  2:  The  second  order  MUSCL  scheme:  density,  N  =  220,  time  t  =  0.3.  (+)  -  numerical 
solutions,  solid  line  -  exact  solutions. 
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Figure  3:  The  third  order  ENO  scheme:  density,  N  =  220,  time  t  =  0.3.  (+)  -  numerical 
solutions,  solid  line  -  exact  solutions. 
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